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NUMERICAL ANALYSIS OF FLOW AND PRESSURE FIELDS IN 


AN IDEALIZED SPIRAL- GROOVED PUMPING SEAL 
by John Zuk and Harold E. Renkel 
Lewis Research Center 

SUMMARY 

A computer program is presented for finding the flow and pressure fields in a 
spiral- grooved pumping seal model for the limiting case of zero clearance. The govern- 
ing nonlinear partial differential equations are solved numerically using the method of 
finite differences and over-relaxation. The program is written in FORTRAN IV and is 
completely described including the program listing, flow charts, and sample problem. 
The computer program calculates the net volume flow rate, velocity profiles, and pres- 
sure distributions for specific axial pressure gradients, Re5molds numbers, and aspect 
ratios. Other biharmonic problems may be solved using this program. 


INTRODUCTION 

In a companion paper (ref. 1) an analysis is given for the flow and pressure fields 
in a spiral- groove pumping seal model for the limiting case of zero clearance. A spiral- 
groove face seal (ref. 1) is a member of a general class of pressure generation devices 
that are characterized by two surfaces moving relative to each other with very small 
film thicknesses and with one or both surfaces grooved. Several geometric forms are 
found; for example, the cylindrical form (viscoseal), the herringbone groove bearing, 
and the conical and spherical bearing forms. Variations and combinations of these are 
also found. The numerical solutions of the exact governing equations using the numeri- 
cal analysis and computer program presented herein are compared with classical models 
of the groove axial flow which neglect the coupling of the groove cross flow. The solu- 
tions presented in reference 1 included the following results: The cross flow shifts the 
pumping flow toward the land leading edge. Conditions under which the classical models 
give good approximations for the relation between axial pressure gradient and net volume 
flow are shown to depend on the Reynolds number and aspect ratio. The groove cross- 



section static pressure is nearly constant except near the moving surface region. A low 
pressure region suggests the possibility of degassing and cavitation; a high-pressure 
region near the land leading edge results in a lift force acting on the moving surface. 

The objective of this report is to present a method of solution and to present a com- 
puter program for numerical solutions for the flow and pressure fields in a spiral- 
grooved pumping seal model whose analysis is given in reference 1. Also, other physi- 
cal problems that can be solved by the computer program are discussed. The computer 
program is written in FORTRAN IV for the Lewis Research Center IBM 709411/7044 
direct-couple system. 


NUMERICAL ANALYSIS 


Seal Model and Equations 


As described in reference 1, the spiral-groove pumping seal model is a stationary 
rectangular cross-section groove with a wall (upper plate) moving at an oblique angle to 
the groove edges as illustrated in figure 1. In addition, a pressure gradient is imposed 
in the groove axial direction. A rectilinear Cartesian coordinate system is used. The 
flow is fully developed in the groove axial direction (z* -direction); that is, end effects 
are neglected in the z*-direction. The flow studied is for a homogeneous, incompres- 
sible Newtonian fluid under steady laminar flow conditions. 

In reference 1 the flow field variables and resulting equations were nondimension- 
alized. (All symbols are defined in appendix A including dimensionless scaling values. ) 

In order to facilitate numerical analysis, the flow field equations across the groove 
(x*-y*-plane) are expressed in terms of a cross flow stream function and the 

vorticity ?*(x*, y*). The derivatives of the stream function are related to the groove 
cross flow velocity components u* and v* such that 




A 


u’^ = 


9y* 


V* = , MU 

3x* 




( 1 ) 


In figure 2, the groove cross flow plane and the stream function direction are shown. 
Since the flow is fully developed in the z*-direction (groove axial direction), 9w*/0z* = 0 ; 
thus, the use of the stream function automatically satisfies the dimensionless incompres- 
sible continuity equation 


2 


( 2 ) 


9u* ^ dv* 
dx* By* 


^cot^a^ = 0 
dz* 


The component of vorticity in the groove axial direction (z*-direction) is 


- X 9v* - 1 
dx* X dy* 



1 d^ip*\ 

^ 0y*V 


(3) 


Two important dimensionless parameters were found from nondimensionalizing the gov- 
erning equations: 


^ bUsincv 

V 


and 


X = d/b 


Groove Cross Flow Plane (x*-y* plane) 

For the stated restrictions the appropriate flow field equation in the groove cross 
flow plane is the two-dimensional vorticity transport equation (see ref. 1) which reduces 
to 


d\* ^ 1 d\* 

3 *2 -2 - *2 
9x* X 9y* 


Re 


dx}/* 9^* dx}/* 

dy* 9x* 9x* 



(4) 


Groove Axial Flow Direction (z*- direction) 

For fully developed flow in the z*-direction, the Navier-Stokes equation is 


dx}/* 9w* _ 9i^ 9^ ^ _ 9P* 1 ^ 9^w* ^ X ^ 9^w* ' 

.2 .2 - * 2 , 

X dy* 


3y* 3x* 0x* 3y* 


+ 

0z* Re^ 


(5) 


.9x*' 


where 
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= Cl = Constant 

dz* ^ 


Hence, 


P* = Cj^z* + C2(x*,y*) (6) 

The boundary conditions will now be stated: For convenience the stream function is 
chosen to be zero on the walls, hence, 

0) = 0 t//*(x*, 1) = 0 y*) =0 ^(1, y*) = 0 (7) 

The fluid velocity no -slip and impermeability condition on the walls expressed in stream 
function form is 

^(x*,0) = 0 i?^(x*,l)-l ^{0,y*) = Q l^(l,y*) = 0 (8) 

dy* dy* 9x* 9x* 

And the no -slip condition for the groove axial direction velocity is 

w*(x*,0) = 0 w*(x*, 1) = -1 w*(0, y*) = 0 w*(l, y*) = 0 (9) 

Once the flow field is found in the x*-y* plane, the static pressure field can be 
found from the dimensionless Navier Stokes equations modified in the following way. 
using the dimensionless stream function and vorticity: 

For Re > 1 (ref. 1): 


9P* 

9x* 


9P* 

9y* 


9y* 9y* 9x* 9x* 


9y*2 


^ d^}p* j 

11 

9?* 

X^* 

dx*^ / 

Re X 

9y* 

dx* 


1 +J-X 

9C* 

xc* 

0x* dy* j 

* Re 

9x* 

'-b 

9y* 


( 10 ) 


( 11 ) 


For 0 < Re < 1 (as stated in ref. 2, the pressure has to be rescaled for small val- 
ues of Reynolds number): 

1 9P*’ I ^ ^2 ^ I _ J_ 1 ^ 

9y* 9y* 9x* 


9x* a *2 
9x* 


Re X 9y* 
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Re 9x* 


9y* 9x* 


9x^ 


( 12 ) 



- x<* 


(13) 


JL = - / 3i//* ^ d^x{/* \ ^ J_^ 3g* 

Re 9y* I 3y* gy*2 9x* 9x* 9y* y Re 9x* 


d\f/* 

9y* 


where P*’ = ReP*. 

In reference 2 results of the integration of the pressure distribution on the moving 
wall surface which yields a net lift force are shown. This net lift force per axial length 
is found from 


-p* 

Axial length 



P* dx* 


(14) 


In reference 2 it was further stated that if degassing occurred at the trailing edge (this 
region would be at ambient pressure) then a net lift force would occur only along the 
leading edge interval. The leading edge interval is defined as the range of x* values 
from the point (where P* is always greater than 0) to the point x* = 1, 


Leading- edge force 
Axial length 


r 

2k 


P* dx* 


P*>0 


(15) 


Outline of Solution 

Due to the geometrical configuration and the nonlinearity of the flow field equa- 
tion (4), analytical solutions are extremely difficult to obtain; however, equations (3), 
(4), (10), and (11) can be solved numerically. 

The basic nondimensional flow field equations (3) and (4) are solved for the stream 
function xj/* and vorticity ?* distributions using finite difference techniques and suc- 
cessive overrelaxation, similar to that used by Lieberstein (ref. 3). Once \p* and ^* 
are known, the normalized pressure field is calculated from equations (10) and (11), 
using a finite difference scheme suggested by Burggraf (ref. 4). From the results of 
equations (3) and (4) the dimensionless z* -directional flow w* and net volume flow Q* 
can be calculated for specified values of the constant groove axial pressure gradient 
9P*/9z*. Equation (5) is solved by the method of finite differences for the w*-field and 

Q* is calculated from 
z 

Q* = w* dy* dx* (16) 

z ^0 *'0 
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Finite Difference Method 


A grid of mesh points (i, j) is constructed over the positive x*-y* plane (fig. 3) 
with i increasing for decreasing y* values and j increasing for increasing x* val- 
ues. With this mesh, equations (3) to (5) can be developed into appropriate difference 
forms for solution on a digital computer. Central differencing techniques were used in 
the equations whenever possible; however, forward or backward differences were some- 
times used especially at the walls. The FORTRAN IV computer programs are described 
in appendix C. The first program for the numerical solution of an idealized spiral- 
grooved pumping seal model calculates the vorticity and stream function, the u*- and 
V*- velocity profiles, the normalized pressure field, and the net lift forces. The second 
program calculates the w*- velocity profile and dimensionless net volume flow rate 
along the groove axis. 

The initial distribution for the stream function is assumed as a linear function with 
i//*' = 0 at the walls and xp* = -0. 1 at the center of the rectangular groove. For the 
case X = 1, square groove, the stream function contours are squares of constant value. 
Various other distributions including ip* equal to a constant for all interior points were 
examined but were found to be less efficient in that more iterations were required to ob- 
tain convergence. The vorticity initial distribution is calculated at all interior points 
using equation (3) which in difference form becomes 




Ax* 





(17) 


At the boundaries the initial vorticity values are calculated from equations (B3) and (B4). 
In finite difference notation these equations in dimensionless form, become 
(1) Lower stationary wall 


y* 

^ imax, i 


XAy*‘ 


(’^imax, i " ’^max- 1, j) 


(18) 


(2) Left stationary wall 




* 

i, 1 


2X 

r 

Ax*' 


k 1 - 2 ) 


(19) 
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(3) Right stationary wall 


^i, jmax 


2X 

Ax*^ 



jmax-l) 


(4) Upper moving wall 




* 

1, j 


2 


(^i_j - - iy’) 


( 20 ) 


( 21 ) 


The solutions to equations (3) and (4) are calculated in an iterative routine in which 
the stream function or vorticity field is scanned once before entering the other field. In 
the literature, other authors (e. g. , ref, 5), who have solved similar boundary value 
problems, have scanned each field a various number of times (from two to 50) before en- 
tering the other field. The present authors feel that this only overcorrects the values in 
that particular field, especially when a relaxation factor is used in the calculations. In 
the iterative process equation (4), which in difference form is 


Ax* 


H 1. i ■ " 'f-i. i) “ 


Re 


4Ax*Ay* 






(22) 


serves as the basis for the vorticity calculations at the interior points; equations (18) 
to (21) are used for the vorticity boundary values. The values of the stream function at 
the interior mesh points in the iterative process are calculated from equation (17), while 
the boundary values are fixed at p* = 0. 

The successive overrelaxation technique used in this analysis is based on a paper 
by Lieber stein (ref. 3). The basic iterative relaxation equation is 


, n+1 


n 


= 


y2»---.yk) 

f'(yi,y2>-*-.yk) 


(23) 


where y? and yf^^ are either pf ^ or . at the n or n+1 iteration; 
f(yi» y 2 > • • • yk) either equation (17) or (22); f’(y^, y 2 > • ■ ■ , yj^) is the combined coeffi- 
cients of pf . from equation (17) or the combined coefficients of . from equa- 
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tion (22); and w is the relaxation factor. In evaluating y 2 > , yj^) and 

f'(yj> 5 ^ 2 ? • • • > Yif) most recent available values for the y*s are used. Substituting the 
appropriate terms into equation (23) yields the following equations for the vorticity and 
stream function as coded in the computer program: 


«r,T‘ ' “ ■ 


MJ+I MJ-I , M-I,i -"M-HJ 


Re 


Ax*" 


^2 . *2 
X Ay* 


4Ax* Ay* 






J 


^+n+l ^ (1 _ ^ 

^1,3 ' '^^1,3 A 




(24) 


X Ay=* 


/i-l,] ^1+1,3, 


+ If*'* 

^i,3 


n +1 


(25) 


where 


K 


1 ~ 





X^ Ay* 7 


Values of w between 0 and 2 are used as relaxation factors, but the optimum value 
for most rapid convergence depends on the aspect ratio and mesh size. For co = 1, this 
iterative process reduces to the Liebmann iterated forms as used by Mills (ref. 5). 

Once the stream function field has been found the u* and v* velocity components in 
the X*- and y*-directions, respectively, can be calculated from equation (1). 
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The stream function boundary conditions (eqs. (7) to (9)) specified in terns of veloc- 
ity boundary conditions are 


u*(0, y*) = u*(l, y*) = u*(x*, 0) = 0, u*(x*, 1) = 1 


v*(0,y*) = v*(l, y*) = v*(x*, 0) = v*(x*, 1) = 0 


(26) 


In central difference notation these equations become 




2 Ay* 




and 




V* . 
1,1 


iL* . ^ - xL* . , 
2 Ax* 


(27) 


These equations are used in this form to calculate the velocities at all interior mesh 
points. 

The static pressure field is obtained by a point-to-point integration of equations (10) 
and (11) over a single mesh width. At each mesh point the values obtained from the two 
equations are averaged to find one pressure for each point. Substituting equation (1) 
into equations (10) and (11) and rearranging terms yields the following forms from which 
the static pressure is calculated: 


= . J_ 1 + X5*v. - u» iV 

dx* Re X ay* 9y* 3x* g^^2 


iZ! . A X ^ - u. ^ . x^v ^V - 

ay* Re ax* gy^2 3x* dy* 


(28) 


(29) 


For the case 0 Re < 1, the corresponding static pressure (P*') equations are obtained 
by multiplying the right side of equations (28) and (29) by Re. 

In the finite-difference representation of these equations the pressure partial deriva- 
tives are approximated using forward or backward differences and the values of the pres- 
sure at adjacent mesh points. Each of the terms on the right side is calculated at the 
same pair of adjacent points and then corresponding terms are averaged to give one ap- 
proximation for each term. For example, in equation (28), the pressure term is approxi- 
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mated using points (i, j) and (i, j+1) on the grid. Each term on the right side is then eval- 
uated at (i, j) and at (i, j+1); the value of a term from (i, j) is averaged with the value 
from (i, j+1) to give one approximation in the finite -difference representation. 

To adjust for the discontinuity in the vorticity at the junctions of the moving wall 
with the stationary walls the following scheme was adopted for these two points: 

(1) If the term being defined is a partial derivative, the value of the vorticity is 
chosen in the same direction as the derivative; for example, the values for 9^*/9x* are 
calculated from equation (21), and the values for d^*/dy* are calculated from equations 
(19) or (20). 

(2) If the term being defined is the vorticity itself, the value chosen is dependent on 
the direction of the pressure derivative term; for example, the value of the vorticity 
used in equation (28) is calculated from equation (21), and the value in equation (29) is 
calculated from either equation (19) or (20). The finite-difference equations representing 
equations (28) and (29) at the interior mesh-points, along the walls, and at the corners 
are thus given by 


(1) Interior points 


Ax* ReX 2 I 2 Ay* 2 Ay* j 2 \ ^ ^ i>J-l 1.3-1/ 


u* 

1,3 


■^tj-l 


*^ 1-1, j+1 ~ ^1 +1, j- 1 ’^x+1, j- 1 ■ "^i+l, j+1 

4 Ax* Ay* 


4 Ax* Ay* 


V* . 

1,3 


fxl/* . . - 2\L* .+;//*. , 
V-3+1 ^1,3 ^i,]-l 


Ax*' 




'\p* . - 2i//* . 1+1/.*. o 
V3 ^1,3-1 1,3-2 


Ax* 


(30) 
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^ ^j~^i+l,3 _ X ^1, j-1 ^i+l,i+l~^i+l, j-l \ 


Ay* 


Re 2 \ 2 Ax* 


2 Ax* 




/\L* , . - 2\pf . + d/* . . 
u* . I 1-1,3 1,3 1+1,] 


1,3 


Ay=* 


-I- 1 

1+1,3 


l//* . - 2xb^ - . + rt . 

'*' 1 , ] *^1+1, ] '^x+2 , ] 


Ay* 


, , n- 1, j*i -n-iA-i*n.i.i-i-n*i.i.i 

3 V 4 Ax* Ay^ 


+ ^ . 

1+1, j 


'n, ni: n,i-i^ 

4 Ax* Ay* 


(31) 


(2) Boundary points 

In equations (32) through (47) which follow, the zero-valued terms due to the boundary 
conditions of the spiral- grooved pumping seal have been omitted, 

(a) Stationary wall (x* = 0, 0 < y* < 1) 


P? o - Pf 


2 ~ "i, 1 _ 1 1 pt- 1, 1 ' ^i+1, 1 ^i- 1, 2 ~ ^i+1, 2 ) X _ "i,2 

9 i, 2 i, 2 « 


uf 


Ax* 


ReX 2 \ 2 Ay* 


2 Ay* 




ti,3-n-i,2^n^i,2-n^i,3 \ ^ ,2 n,2 n, 4- 3^ n, 2 


2 Ax* Ay* 


Ax* 


(32) 
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( 33 ) 


P* 


X Ml, 2 -^. 1 , 1 


Ay=* 


Re 2 \ Ax* 


Ax* 


(b) Stationary wall (x* = 1 , 0 < y* ■< 1 ) 


p* _p* ^ K* - r* if* +1°* 1 

i, jmax i, jmax-1 _ _ 1 i pi-l, jmax ^i+1, jmax ^ ^i-1, jmax-1 ^i+1, jmax-1 

Ax* ReX 2 I 2 Ay* 2 Ay* 

+ ^f* V* jmax -1 

n^i, jmax-1 i, jmax-1 « 


'^i-1, jmax- 1 ~ 1, jmax-2 ’^i+l, jmax -2 ’^i+ l, jmax-1 

2 Ax* Ay* 


o v?‘ . 1 /li'? . ^ . n+xpf . o 

^ ^2 1, ]max-l ^1, ]max-l ^1,3 max - 2 ]max -3 


( 34 ) 


Ax*' 


p» _ p* /if* _i"* i°* _if* 

i, jmax i+1, jmax _ 1_ ^i, jmax jm ax- 1 ^ ^i+1, jmax ^i+1, jmax-1 

Av* Re 2 \ Ax* Ax* 


( 35 ) 


(c) Stationary wall (y* = 0 , 0 < x* < 1 ) 


p* _ p* 

imax, j imax, j- 1 _ _ 1 _1 


Ax* 


ReX 2 


X 


^imax-1, j ~ ^imax, j ^imax- 1, j- 1 ~ ^imax, j- 1 


Ay=* 


Ay=t 


( 36 ) 
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p* p+ -If* 

i, max- 1, i imax, j X i pimax, j+1 ^Imax, j-1 ^ ^imax-l,j+l ^imax-l,j-l 


Ay* 


Re 2 


'2 Ax* 


2 Ax* 


hf* u* *^imax- 1, j 

2 ^imax- 1, j imax- 1, j 2 


^ I ’^imax-3, j ~ ^’^imax-2, j ^imax-1, j | + ;y2 ^imax- 1, j 

Ay*2 / 2 


’^imax-2, j+1 " ’^im ax-2, j-1 ’^ imax-l,]-! ~ ’^imax-1, j+1 

2 Ax* Ay* 


(37) 


(d) Moving wall (y* = 1, 0 < x* < 1) 


Ax* ReX 2 \ Ay* Ay* 


2 \ 2 Ax* Ay* 


. n,rn.i-2*n.i-2-n,i 

2 Ax* Ay* 


(38) 
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Ay* Re 2 \ 2 Ax* 2 Ax* J 2 I J J 3 


1 ... Aii 

2 y Ay»2 y y 




,.2Vll Ki-.l-»ii-l"»ji,l- l- 'fS.i+l 

2 \ 2 Ax* Ay* 


(e) Corner point (x* = 0, y* = 1) 


^ 1 , 2 ~ ^ 1 , 1 _ 1 1 , ^ 1,2 ~ ^^2 


ReX 2\ Ay* 


1 m, 2 ~ 1 ^^2, 1 ~ *^2, 2 

2 1 Ax* Ay* 


n,3-n,2^H,2~n 

Ax* Ay* 




Re 2 V Ax* 


(f) Corner point (x* = 1, y* = 1) 
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imax~ ^l,jmax-l ^ 1 jmax ^ ^1, jmax-1 ~ jmax-1 


Ax* 


ReX 2 \ Ay=* 


Ay-1 


_ lj ’^1, jmax ^1, jmax-1 ''' jmax-1 ~ jmax 

Ax* Ay* 


^ *^1, jmax- 1 “ jmax-2 ^^2, jmax-2 ~ ^2, jmax-1 

Ax* Ay* 


(42) 


^1, jmax~ ^^,jmax _ 1 p 1, jmax ~ ^1, jmax-l ^ jmax ~ jmax-1 

Ay* Re 2 \ Ax* Ax* 


(43) 


(g) Corner point (x* = 0, y* = 0) 


^imax, 2 ~ ^imax, 1 _ _ 1 1 / ^imax-l, 1 ~ ^imax, 1 ^ ^imax- 1, 2 ~ ^imax, 2 

Ax* ReX 2\ Ay* Ay* 


(44) 


p* _ p* - ^* If* - if* 

imax-l, 1 imax, 1 _ X 1 ( ^imax, 2 ^imax, 1 ^imax-1, 2 ^imax-1, 1 


Ay=* 


Re 2 \ Ax* 


Ax* 


(45) 


(h) Corner point (x* = 1, y* = 0) 

P* - P* ^ It* - C* 

imax, imax imax, jmax- 1 _ _ 1 1 pimax-1, jmax ^imax, jmax 

Ax* ReX 2 \ Ay* 


r* _ if + 

^ imax- 1, jmax- 1 ^ imax, jmax-1 

Ay* 


( 46 ) 
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p? . . - p* 

imax- 1. jmax imax, ]max 

Ay* 


\ 1 
Re 2 


If* 

^imax, jmax 


If* 

^imax, jmax-1 


Ax* 


+ 


jf * 

^iinax- 1, jmax 


4max- 1, jmax-1 


Ax* 


(47) 


The initial distribution across the entire static pressure field is constant at P* = 1. 

The program then iterates on the pressure field from top-to- bottom and left-to-right un- 
til convergence is achieved. The normalized static pressure field is calculated by sub- 
tracting from all the pressures in the field the value at one reference point specified as 
a program input. It is felt that this normalization method presents more useful results 
rather than the pressure ratios presented in references 1 and 2. 

Both integrations in equations (14) and (15) in finding the net lift force per axial 
length and leading edge force were performed numerically using Simpson's rule. 

In program two, the w*-veIocity profile is calculated based on the stream function 
distribution calculated in program one. The stream fimction data are read in to pro- 
gram two, and the w*-field is initially equal to one, except for the boimdary conditions 
presented with equation (5). If for a certain set of parameters. Re and A, the w*-field 
is to be calculated for more than one value of 9P*/3z*, the initial value of the w*-field 
for the succeeding 9P*/9z* value is set equal to the converged value from the preceed- 
ing 9P*/9z*. Lieberstein’s (ref. 3) successive over relaxation technique is again used 
in the solution of equation (5) with the y.'s of equation (23) being equal to the latest 
available values of wf j. The f(y^, y 2 , . . . , yj^) is evaluated from the following equation 
which is the finite difference form of equation (5): 


2 Ay* l\ 2 Ax* 



j+1 ' j- 1\ ri 1, j ~ j \ ^ 9P* 


2 Ax* 


2 Ay=* 


9z* 


J_ 

Re 


Ax*^ 


1 


^i,j ^^i+l ,i 

Ay*^ 


= 0 


(48) 


And f'(yi,y<>, . . . ,y ) is the combined coefficients of wf . from the preceding equation, 
i z n 1, J 
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Substituting the appropriate expressions for f and f into equation (23) yields the 
z*-directional flow equation as coded in program two: 




n+1 


= (1 - a))w?” + Re 


’^ i-1, j ' ’^i+l, 3 
2 Ay* 


rw? 




2 Ax* 


n+1 





(49) 


where is defined in connection with (24) and (25). The range of the relaxation factor 
w is between 0 and 2. 0 with an optimum value of approximately 1.3 for X = 1 and 
29 mesh points in the x*- and y*-directions. Other combinations of X and grid size 
have different optimum relaxation factors. 

After the w*-field has iterated to convergence within the prescribed error condi- 
tion, the net volume flow is calculated from equation (16) using the method of mechanical 
cubature. A requirement of this method, which is based on Simpson's Rule of integra- 
tion for two dimensions, is that there be an odd number of mesh points in both the x*- 
and y*-directions so that the proper weight factors will be applied in the cubature 
scheme. The double integral in equation (16) is thus approximated by the double sum- 
mation as given by 


imax- 1 




_ Ax* Ay* 
9 


i: 


i=2,4,6,... j=2,4,6,... 


jmax-l 

Z I'n, j 


■"'"r-i.i-i 


(50) 
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Convergence Remarks 


A solution was considered to have converged when the following criteria were met: 

1. The relative change in the values of the stream function, pressure distribution, 
or axial velocity distribution between two successive iterations was less than a pre- 
scribed maximum for all mesh points; for example, 

m+1 m 

^ — < Prescribed maximum 

^n+l 

2. The values of the stream function and vorticity at any mesh point changed less 
than 1 percent when the number of mesh points was doubled in either direction. 

In addition to the convergence criteria, the effect of round-off error on the results 
of reference 1 was checked by doubling the precision of the computing machine calcula- 
tions. 


Computer Program Formulation 

As mentioned in the previous section, the equations in finite difference form are 
solved on a high-speed digital computer. A complete description of the program is pre- 
sented in appendix C. The program listing is given in appendix D. Figures 4 to 12 pre- 
sent the computer program flow charts. Flow charts of the main program, stream func- 
tion and vorticity, u* - and V*- velocities, pressure field, w*-velocities, and net volume 
flow Q* calculations are shown. A sample problem for the case where X = 1, 

Re = 100, and (9P*)/(0z*) = -0. 115 with its input and output is given in appendix E. 
(Plots of this case are in ref. 1. ) 


Program Use for Other Physical Problems 

The computer program finds the stream function and vorticity in the groove cross- 
flow plane by solving equations (3) and (4), which are the following: 




(3) 


= Re I ^ 

0x*^' X^' 0y*^ 


(4) 
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These equations can be combined and placed in the following form (The aspect ratio 
X can be eliminated by redefining the independent variables): 

, 4 ^. , ( 51 ) 

0(x*,y*) 

The computer program can be used to solve many problems in mathematical physics that 
appear in this form or are reducible to this form. For example, the biharmonic equa- 
tion (This is the creeping flow case discussed in ref. 1) V i/-* = 0 

Laplace's equation 


= 0 

Poisson's equation 

Of course, the proper transformation of the variables must be made to the form 
solved in the computer program. The proper boundary conditions must be specified. 

The coding in the FORTRAN IV computer program contains all the boundary terms 
including the zero-valued terms omitted in equations (32) through (47) which were spe- 
cifically written for the spiral-grooved pumping seal. The configuration must be rec- 
tangular. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, October 29, 1970, 

126-15. 
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APPENDIX A 


SYMBOLS 


b 

C 

d 

F 

F* 

n 

P 

p*t 

p* 

3P* 

az* 

Q 

Q* 

Re 

U 

u 

u* 

V 

V* 

w 

w* 

X 

X* 

Ax* 

y 

y* 


groove width 
numerical coefficient 
groove depth 

net moving surface lift force 

dimensionless lift force, (F/bpU sin^ a) X (total axial length) 
numerical coefficient 
iteration number 
static pressure 

dimensionless pressure, Pb/pU sin a (for 0 < Re < 1) 

2 2 

dimensionless pressure, P/pU sin a (for Re > 1) 
dimensionless groove axial pressure gradient, constant 

net volume flow rate, groove axial direction 
dimensionless net volume flow rate, Q /(dbU cos a) 

Reynolds number, (bU sin a)/v 
moving wall velocity 
velocity in x-direction 
dimensionless velocity, u/(U sin a) 
velocity in y- direction 
dimensionless velocity, b/djy/(U sin 
velocity in z- direction 
dimensionlesss velocity, w/(U cos a) 
cross groove coordinate 
dimensionless coordinate = x/b 
mesh width in x*-direction 
groove depth direction coordinate 
dimensionless coordinate = y/d 
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Ay* 


z 

z* 

a 

X 

M 

u 

0) 


a(x*, y*) 
Subscripts: 
i 

imax 

3 

jmax 

0 


mesh width in y*-direction 

groove axial coordinate 

dimensionless coordinate, z/(b tan a) 

angle between moving wall direction and groove axis 

aspect ratio, equal to groove depth to groove width, d/b 

absolute viscosity of fluid 

kinematic viscosity of fluid 

stream function, defined in x-y plane 

dimensionless stream function, sin a) 

vorticity component in z-direction 

dimensionless vorticity, ?b/(U sin a) 

relaxation factor 


Laplacian operator, { A 


+i 


ax*^ ^ ay*^ 


biharmonic operator, ( X 


a^ +2 


ax*^ ^ ax*^ ay*^ a ay*'^ 


Jacobian operator, 


dt}/* av^i//* d\j/* dv^xi/* 

ax* ay* ay* ax* 


mesh point in y*- direction 

maximum value of mesh point in y* -direction (y* = 0) 
mesh point in x*-direction 

maximum value of mesh point in x* -direction (x* = 1) 
boundary reference point 


APPENDIX B 


FIRST ORDER VORTICITY BOUNDARY VALUE APPROXIMATION 

The boundary conditions for the stream function are known on the boundary because 
the normal and tangential velocities must satisfy the impermeability and no-sIip condi- 
tions at the wall. The vorticity value on the boundary must be calculated and will vary 
from iteration to iteration until a specified convergence criterion is satisfied. The vor- 
ticity at the boundary can be found by the following first order approximation. 

Consider the boundary along the x*-axis as shown in figure 13. The stream function 
expanded in a Taylor series about point 0 is 



Using equation (3) which relates the stream and vorticity functions yields 



Since ip* = ip*(x*) on the boundary shown in figure 13, 



Substituting equation (B2) into the truncated Taylor series expansion (Bl) results in 




* - 
0 ■ 


2 

X(Ay*)^ 


(Ay*) 



(B3) 


In like manner if the wall were along the y*-axis, the wall vorticity equation is 




* - 
0 " 


2X 

(Ax*)^ 


xl/*-xP* + (Ax*) 



(B4) 
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The resulting vorticity boundary conditions are found by substituting the boundary 
conditions equations (7) and (8) into equations (B3) and (B4) and are for stationary walls: 

^ (horizontal) 

X(Ay.)2 

-2X\p* 

^6 (vertical) 

(Ax*)^ 


and for a moving wall 



2(-;//* + Ay*) 
\(Ay*)^ 


where is the stream function value in the flow field which is Ax* or Ay* away 
from the boundary whose stream function has a value (See fig. 13, ) 
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APPENDIX C 


DESCRIPTION OF COMPUTER PROGRAMS 
Program ISGPSM 

The computer program for the numerical solution of an idealized spiral grooved 
pumping seal model (ISGPSM) consists of the MAIN program and the following sub- 
routines: 

SVCALC enters the input data and calculates the stream function and vorticity distribu- 
tions and many of the constant terms used throughout the program. 

INDIST generates the stream function initial distribution. 

ZETAF generates the vorticity initial distribution. 

VORTWL calculates the values of the vorticity along the moving and stationary boundary 
walls. 

UVFUNC computes the u* and v* velocity profiles in the x* and y* directions, 
respectively. 

PRESS calculates the static pressure field based on the stream function and vorticity 
distributions and then the normalized pressure field with reference to a predeter- 
mined point. 

BCDUMP punches data in sequentially numbered absolute binary cards with a maximum 
of 22 words per card. 

Program ISGPSM is structured to make use of the overlay feature of the computing 
machine loading system (IBLDR) which saves in auxiliary storage those sections of the 
program currently not being executed. The use of this feature thus provides for a maxi- 
mum of 65 mesh points in both the x*- and y*- directions. Figure 14 is a chart of the 
overlay structure of the entire ISGPSM program. 

The input to program ISGPSM is a comparatively few number of variables which are 
entered on one data card. The following is a list of the variables and the definition of 
each: 

LAMBDA aspect ratio 
RE Reynolds number 

OMEGA relaxation factor [0, 2. 0] 
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PCT 

NOPTY 

NOPTX 

NOLESTY 

NOLINX 


maximum allowable relative change between two successive iterations of the 
stream function or static pressure field calculations at any mesh point 

number of mesh points in y*-direction (maximum number = 65) 

number of mesh points in x*-direction (maximum number = 65) 

number of lines of output data in y*-direction (maximum number = NOPTY) 

number of columns of output data in x*- direction (maximum number = 15) 


NOLINY and NOLINX should be chosen such that (NOPTY - l)/(NOLINY - 1) = and 
(NOPTX - l)/(NOLINX - 1) = K„ where K and K are integers. The program will 

X X y 

then print lines 1,1 + K^, 1 + 2Ky, . . . , NOPTY and columns 1,1+ K^, 1 + 2K^, . . . , 
NOPTX. 


IREF mesh point number in the y*-direction for predetermined reference point in 

the normalized pressure field calculations 


JREF mesh point number in the x*-direction for the predetermined reference point 

in the normalized pressure field calculations. 


NOW a control word for the w*-velocity calculations; if NOW equals any nonzero 

integer, the converged values of the stream function at all mesh points 
will be punched into cards for use in w*- velocity profile program, but if 
NOW equals zero, no card punching will occur. 

NOPLOT a control word for the normalized static pressure field plots. If NOPLOT 
equals any nonzero integer, the normalized static pressure field is 
punched using program BCDUMP for use with Canright and Swigert (ref. 6) 
three dimensional plotting program, but if NOPLOT equals zero, no card 
punching will occur. 


PROBNO a control word for the tjqje of problem to be solved. If PROBNO ^ 0, the 
program will solve Laplace's equation (£* = 0) or Poisson's equation 
(f* = f(x, y)). If PROBNO = 0 and Re = 0, a solution to the biharmonic 
equation is obtained. 


FLUFLG a control word for additional fluid flow calculations. If FLUFLG = 0, the 
u— and V— velocity profiles and pressure fields will be calculated. 


The format for this card is (4F8. 0, 1013). 


16 

17 

24 

25 

32 

35j 

38 

41 

44 

47 

50 

53| 

56 

59 

62 

XXX. 


x.xx 


.XXX 

XX 

XX 

XX 

XX 

XX 

XX 

X 

X 

X 

X 
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The printed output from the ISGPSM program includes the aspect ratio, Reynolds 
number, relaxation factor, maximum allowable relative change, number of iterations to 
convergence of the stream function, and the number of increments in the x* and y* 
directions. Also printed are paragraphs of data NOLINY lines long and NOLINX columns 
wide of the stream function and vorticity, the u* and v* profiles, and the static and 
normalized pressure fields with the number of iterations to convergence of the static 
pressure field. - Following the normalized pressure field are the leading edge region and 
whole surface net forces. 


Program WFIELD 

Program WFIELD is a FORTRAN IV computer code for calculating the z*-direction 
velocity profile w* and the dimensionless net volume flow rate Q* along the groove 
axis. LAMBDA, RE, OMEGA, PCT, NOPTY, NOPTX, NOLINY, and NOLINX are the 
first input variables to WFIELD and are entered on one data card. The definition of these 
variables is the same as for the first eight variables in the ISGPSM program and the for- 
mat for this card is (4F6. 0, 413) 


1 6 

7 

12 

13 

18 

19 

24 

27 

30 

33 

36 

37 

72 

73 

80 

X.XXX 


XXX. 


X.XX 


.XXX 

XX 

XX 

XX 

XX 






The second set of input cards is the values of the stream function (PSI) at all mesh 
points. These cards are punched in the format (9F8. 2) by the ISGPSM program when 
the variable NOW is not equal to zero. The number of cards required is dependent on 
the grid size with a maximum of nine data words per card. The dimensionless groove 
axial pressure gradient 9P*/8x*, which in the program is denoted by DPDZ, is the last 
input variable. For a specific geometry and stream function distribution, the w*- 
velocity profile and net volume flow may be calculated for several 9P*/9z* values; 
however, each value must be entered on a separate data card in the format (F8. 0). 


1 8 

9 

72 

73 

80 

X.XXX 
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The output from the WFIELD program is the aspect ratio, Reynolds number, relaxa- 
tion factor, maximum relative change, number of iterations to convergence, number of 
increments in the x* and y* directions, and pressure gradient. Associated with a 
particular pressure gradient, the w*- velocity at the specified mesh points and the net 
volume flow along the groove axis are also printed. 


Program GRAPH 


The three-dimensional plots of the normalized pressure field are generated using 
program GRAPH, subroutine PSURF, and the PLOT3D package of subroutines of Canright 
and Swigert (ref. 6). (tt is assumed that the users of program GRAPH have the appro- 
priate Calcomp plotter available although this is not necessary to use programs ISGPSM 
and WFIELD. ) The data to be plotted is read by the computer via the FORTRAN IV pro- 
gram GRAPH. The PLOT3D subroutines then analyze the data to establish the array of 
coordinates of each point to be plotted, scale these values to the size of the plotting 
paper, and set up the figure axes. If the axes are to be rotated to present a nonstandard 
three dimensional projection, the data points and figure axes are transformed to the new 
coordinate system. The information on the data points, figure axes, and figure labels is 
then written on a magnetic tape for further processing by a California Computer Products 
(Calcomp) magnetic tape plotting system. 

The input to program GRAPH is LAMBDA, NOPTY, NOPTX, and 


NOYPLS 

NOXPLS 


number of planes parallel to x*-P* plane or perpendicular to the y*-axis 
number of planes parallel to y*-p* plane or perpendicular to the x"*=-axis 


XSCALF 

YSCALF 

PSCALF 


( X*-, y*- and P*-direction scale factors used to adjust the coordinates of 
each point to match the size of the plotting paper and present the plot in 
the proper perspective (These scale factors are dependent on the amount 
. of hardware on the Calcomp plotter. ) 


These variables are entered on one data card in the format (F6. 0, 413, 3F6. 0). 


1 6 

7 9 

10 

12 

13 

15 

16 

18 

19 

24 

25 

30 

31 36 

37 

72 

73 

80 

X. XXXi 

XX 


XX 




XX 


X.X 


X.X 

X,J^ 
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In addition, the normalized pressure distribution at all mesh points is read from punched 
cards via the BCREAD input routine. These cards are generated in subroutine PRESS of 
the ISGPM and are punched by the BCDUMP routine in column binary format. 


1 
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APPENDIX D 


PROGRAM LISTING 


SIBFTC MAIN LIST, DECK 
C 

COMMENT — NUMERICAL SOLUTIONS OF CONVECTIVE INERTIA EFFECTS FOR AN 
C IDEALIZED, RECTANGULAR, PARALLEL GROOVED PUMP - SEAL MODEL. 

C 

COMMON PSI ( 65,65) ,2ETA( 65,65) ,NOPTX,NOPTY,NPTXM1,NPTYM1 ,DX,DY, 

1 TOOOXS,TOOOYS,LAMBOA,TWODX,TWODY,RE,XINDX,YINOX,PCT, IREF, JREF, 

2 U( 65,65) ,V( 65,65) , NOW, NOPLOT , F LUFLG, PARTS I , LMODXS , LMDYSO, PROBNO 
INTEGER FLUFLG 

1 CALL SVCALC 

IF (FLUFLG ,EQ. 0) CALL PRESS 

GO TO 1 

END 


$IBFTC CALCSV LIST, DECK 
SUBROUTINE SVCALC 
C 

COMMENT — STREAM AND VORTICITY DISTRIBUTIONS. 

C 

COMMON PSI (65,65) ,ZETA( 65,65) , NOPTX , NOP T Y , NPT XMl , NPTYMl , OX , DY , 

1 TOODXS,TOODYS, LAMBDA , TWODX, TWODY,RE,XI NDX , Y I NDX , PCT , IREF, JREF, 

2 U( 65,65) ,V( 65,65) , NOW, NOPLOT, FLUFLG, PARTS I , LMODXS , LMOYSQ, PROBNO 
DIMENSION PSIOUT( 65,65 ) 

REAL LAMBDA, LMODXS, LMOYSQ, LMSOYS 
INTEGER X I NDX, Y I NDX, PROBNO, FLUFLG 
LOGICAL JAIL 
WRITE (6,100) 

100 FORMAT ( 1H18X,38HI0EALIZED SPIRAL GROOVED PUMPING SEAL. /IX) 

75 READ ( 5,3 ) L AMBOA , R£ , OMEG A , PC T , NOPT Y , NOPTX , NOL I NY , NOL I NX, IREF, JREF 
A , NOW, NOPLOT, PROBNO, FLUFLG 
3 FORMATl^FB. 0,1013) 

NPTXM1= NOPTX-1 
NPTYM1= NOPTY-1 
210 XINOX = NPTXM1/(NQLINX-1) 

YINDX = NPTYMl/ (NOLINY-i) 

NPTXPl = NOPTX+1 
NPTYPl = NOPTY + 1 
DX = l./FLOAT(NPTXMl) 

DY = l./FLOAT(NPTYMl) 

DXSQ = OX • DX 
DYSQ = DY • DY 
TWODY = DY + DY 
TWODX = DX + DX 
LMOOXS= LAMBDA/DXSQ 
LMDYSO = LAMBDA*DYSQ 
LMSDYS = LAM8DA»LMDYSQ 
READXY = RE/4./DX/DY 

PSIMLT = LMDYSQ*DXSQ/2./(LMSDYS+0XSQ) 

ZETMLT = LAMBDA*PSIMLT 
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Toonys = 2 ./DysQ 
TOOOXS = 2./0XSQ 

PARTSI =-2.*(LAMBDA/DXSQ + l./LAMBDA/DYSQ) 

PARTZE = PARTSI/LAMBOA 
OMOPSI = OMEGA/PARTSI 
OMOPZE = OMEGA/PARTZE 
TERM * OMOPZE *RE4DXY 
ONEMOM= 1. -OMEGA 

OMMENT — STREAM FUNCTION INITIAL DISTRIBUTION. 

CALL INOIST 
C 

COMMENT — VORTICITY INITIAL DISTRIBUTION. 

C 

CALL ZETAF 
15 ITKONT = 0 
C 

COMMENT — ITERATIVE SOLUTIONS. 

C 

17 ITKONT= ITKONT+1 
JAIL= .FALSE. 

IF (PRORNO .NE. 0) GO TO 21 
CALL VORTWL 
DO 20 J=2»NPTXM1 
DO 20 I=2,NPTYM1 

ZETA(I.J) = ZETAd, J)*ONEMOM - OMOPZE* { ( Z£ TA ( I , J + 1 ) +ZET A ( I , J-1 ) )/ 

1 DXSQ * (ZETA(I+1,J)+ZETA(I-1,J))/LMSDYS) 

20 IF (Rfc .GT. 0.)ZETA(I,J) = ZETA (I , J ) -TERM* d PS I ( I . J + 1 J-PS 1 1 I . J- 1 )) 

1 *(ZETA( I-l,J)-ZETA(I+l,jn - ( PSI (I -1 , J ) -PS I ( I +1 , J )) » ( ZE TA ( I , J+1 ) 

2 -ZETAI I , J-1) ) ) 

21 DO 22 J=2fNPTXMl 
DO 22 I=2,NPTYM1 

PSIF= PSI 1 1 » J)*ONEMOM - OMOPSI * ( LMOOXS* ( PS 1 1 I . J + 1 ) + PS I ( 1 1 J~1 ) ) + 

1 (PSI(I+1,J)+PSI(I-1,J) )/LMDYSQ + ZETAII.J)) 

IF (JAIL) GO TO 22 

23 IF (ABS( (PSIF-PSI ( It J) )/PSIF) .GT. PCT) JAIL = .TRUE. 

22 PSI(ItJ)= PSIF 

IF (JAIL) GO TO 17 

45 WRITE (6,46) LAMBDA , RE , OMEGA , PC T , I TKONT , NPTXMl , NPTYMl 

46 FORMAT (9X.8HLAM8DA =F6. 3 . 5X , 4HRE =£12.4,5X, 1 9HREL AX AT ION FACTOR 
A=F5.2/9X,24HMAXIMUM RELATIVE ERROR =F9. 6 , 5X , 22HNUMBER OF ITERATION 
BS = I5/9X, 37HNUMBER OF INCREMENTS IN X-DIRECTION =I4,5X,16HIN Y-CIR 
CECTION =I4/1H08X,24HSTREAM FUNCTION X 10E+4./1X) 

C 

COMMENT — STREAM FUNCTION AND VORTICITY OUTPUT . 

C 

DO 50 J=1,N0PTX 
DO 50 I=1,N0PTY 

50 PSI0UT(I,J)=PSI(I,J)»l.E+4 

IF (NOW .NE. 0) PUNCH 48 , ( ( PSI OUT ( I , J) , 1=1 , NOPTY) , J= 1 , NOPTX ) 

48 F0RMAT(9F8.2) 

DO 51 1=1, NOPTY, YINDX 

51 WRITE (6,52) ( PS TOUT ( I , J ) , J=1 , NOPTX , XI NDX ) 

52 F0RMAT(9X,15F8.2) 

WRITE (6,64) 

64 FORMAT( 1H18X,10HV0RTICITY,/1X) 

DO 67 1=1, NOPTY, YINDX 

67 WRITE (6,68) ( ZETA ( I , J ) , J = 1 , NOPTX , X I NDX ) 
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68 F0RMAT(9X, 15F8.3) 

IF (FLUFLG .EO. 0) CALL UVFUNC 

RETURN 

END 


$IBFTC DISTI LIST, DECK 
SUBROUTINE INDIST 

COMMON PS I (65,65) ,ZETA(65,65) , NOPTX, NOPT Y, NPTXMl , NPT YM I , DX , DY , 

1 TOODXS,TnODYS,LAMBDA,TWODX,TWODY.RE,XINDX,YINDX,PCT, IREF, JREF, 

2 U( 65,65) ,V(65, 65) , NOW , NOPLOT , F LUFLG , PARTS I , LMODXS , LMDYSQ, PROBN 0 


NPTXPl 


NOPTX+1 

NPTYPl 

= 

NOPTY+l 

MIDPTX 

= 

NPTXPl/2 

DO 8 J 


1, MIDPTX 

PSITRM 

= 

-.2*FL0AT( J-1 )*DX 


MAX = NPTXPl -J 
DO 8 I=1,NQPTY 
PSI ( I , J )=PSITRM 

8 PSI ( I,MAX)=PSITRM 
MIDPTY=NPTYPl/2 
JMIN=0 
JMAX^NPTXPl 

DO 10 I=ltMIDPTY 

PSITRM = ^.2*FL0AT( I^i)*DY 

MAX = NPTYPl - I 

JMIN = JMIN+1 

JMAX = JMAX -1 

IF (JMIN .LE* JMAX) GO TO 9 

JMIN = NPTXPi/2 

JMAX = (NPTXPl + l)/2 

9 DO 10 J=JMIN,JMAX 
PSI ( I, J )=PSITRM 

10 PSI (MAX,J)=PSITRM 
RETURN 
END 


$IBFTC ZETAXY LIST, DECK 
SUBROUTINE ZETAF 

COMMON PSI (65,65) ,ZETA( 65,65) , NOPTX, NOP TY, NPTXMl , NP TYMl , D X , DY , 

1 TOODXS,TOODYS,LAMBDA,TWOOX,TWOOY,RE,XI NOX , Y I NDX , PCT , I R EF , J REF , 

2 U( 65,65) ,V(65,65) , NO W, NOPLOT , FLUFLG , PARTS I , LMODXS , LMDY SQ, PROBNO 
REAL LMODXS, LMDYSQ 

14 DO 16 J=2, NPTXMl 
DO 16 I=2,NPTYM1 

16 ZETA(I,J)=-PARTSI*PSI(I,J)-(LM0DXS*(PSI(I,J+1)+PSI(I,J-*1) ) + 

l(PSI(I+l,J)+PSI(I“ltJ)) /LMDYSQ) 

DO 12 J=l, NOPTX, NPTXMl 
DO 12 I=1,N0PTY,NPTYM1 
12 ZETA( I, J ) =0.0 
RETURN 
END 



$I BFTC VRTWAL LI ST, DECK 
SUBROUTINE VORTWL 
C 

COMMENT — VORTICITY BOUNDARY VALUES. 

C 

COMMON PS I (65,65) , ZETA ( 65 ,65 ) , NOPTX , NOPT Y, NPTXMl , NPTYMl , DX , DY, 

1 TOOOXS,TOOOYS,LAMBDA,TWODX,TWODY.RE,XINDX,YINDX,PCT, IREF, JR6F, 

2 U( 65,65) ,V (65,65 ) ,NOW, NOPLOT, FLUFLGfPARTS I ,LMOOXS,LMOYSQ,PROBNO 
REAL LAMBDA 

DO 10 J=2, NPTXMl 
C 

COMMENT — MOVING WALL. 

C 

ZETA( 1, J) = (P.il < 1, J)-PSI (2, J)-DY)*TOODYS/LAMBDA 

C 

COMMENT — STATIONARY WALLS. 

C 

10 ZETA(NOPTY,J) = ( PSI ( NOPTY,J )-PSI (NPTYMl ,J )) *TOODYS/LAMBDA 
DO 20 1=2, NPTYMl 

ZETAII,!) = (PSI(I,1)--PSI(I,2))*T00DXS»LAMBDA 
20 ZETA( I,NOPTX) = ( PSI ( I , NOPTX) -PSI ( I , NPTXMl )) *TOODXS*LAMBOA 
RETURN 
END 


$IBFTC FUNCUV LIST, DECK 
SUBROUTINE UVFUNC 
C 

COMMENT — U* AND V* VELOCITY DISTRIBUTIONS. 

C 

COMMON PSI(65,65) , ZETA (65,65 )• NOPTX, NOPT Y, NPTXMl, NPTYMl, OX, DY, 

1 TOOOXS,TOODYS,LAM8DA,TWODX,TWOOY,RE,XINOX,YINOXtPCT, IREF, JREF, 

2 U( 65,65) ,V(65,65) ,NOW,NOPLOT,FLUFLG,PARTSI,LMODXS,LMDYSQ,PR06NC 
INTEGER XINDX,YINDX 

101 DO 102 I=l,NOPTY 
U(I,I)=0.0 

V( I,1)=0.0 
U( I ,NOPTX)=0.0 

102 V( I ,N0PTX)=0.0 

DO 104 J=l, NOPTX 
U( 1, J ) = 1.0 

V(1,J) = 0.0 

U(NOPTY,J) = 0.0 

104 V(NOPTY,J) = 0.0 
DO 106 J=2, NPTXMl 
DO 106 1=2, NPTYMl 

U(I,J) = (PSI(I-1,J) - PSI (1+1, J) )/TWODY 
106 V(I,J) = -(PSI ( I, J+1 )-PSI ( I, J-1) )/TWOOX 
WRITE (6,1000) 

1000 FORMAK 1H18X,20HU* VELOCITY PROFILE, /IX) 

DO 1002 I=1,N0PTY,YINDX 
1002 WRITE (6,1004) ( U ( I , J ) , J= 1 , NOPTX , XI NDX ) 

1004 F0RMAT(9X,15F8.4) 

WRITE (6,1007) 

1007 FORMAT( IH18X,20HV* VELOCITY PROFILE, /iX) 

DO 1010 I=1,N0PTY,YINDX 
1010 WRITE (6,1004) ( V ( I , J ) , J= 1 , NOPTX , X I NDX ) 

RETURN 

END 
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SIBFTC PCALC LIST, DECK 
SUBROUTINE PRESS 
C 

COMMENT — PRESSURE FIELD CALCULATION. 

C 

COMMON PSI165,65) ,ZETA( 65,65 ) ,NOPTX, NOP TY,NPTXM1,NPT YMl, OX, DY, 

1 TOODXS, TOODYS,LAMBDA ,TWODX,TWODY,RE,XINDX,YINDX,PCT, IREF, JREF, 

2 U( 65,65) ,V( 65,65 ), NOW, NOPLOT, FLUFLG, PARTS I , LMOOXS, LMOYSQ, PROBNO 
DIMENSION XTERM(65,65 ) , YTERM ( 65 , 65 ) , P ( 6 5 , 6 5 ) , PNORM ( 65 , 65 ) 
EQUIVALENCE ( P , P S I ) , I PNORM , ZETA ) 

INTEGER XINDX,YINDX 
INTEGER PROBNO 
LOGICAL JAIL 
REAL LAMBDA 
REAL LEREG 

FOURDY = TWODY + TWODY 
FOUROX = TWODX + TWODX 
NPTXM2 = NPTXMl-1 
NPTYM2 = NPTYMl - 1 
NPTXM3 = NPTXM2-1 
NPTYM3 = NPTYM2-1 
C 

COMMENT — INTERIOR POINTS. 

C 

CON2 = LAMBDA/2. 

CON3 = TWODX*FOURDY 
CONA = (LAMBDA/DX)**2/2. 

C0N6 = TWODY*DY 

CON7 = LAMBDA**2/CON3 

IF (RE .GT. 1.) GO TO 7 

CON2 = CON2*RE 

CON3 = CON3*RE 

CONA = CONA*RE 

C0N6 = C0N6*RE 

CON7 = C0N7*RE 

CONI = -l./LAMBDA/FOUROY 

CONS = LAMBDA/FOURDX 

GO TO 8 

7 CONI = -l./RE/LAMBDA/FOURDY 
CONS = LAMBDA/RE/FOURDX 

8 DO 10 1=2, NPTYMl 
DO 9 J=3,NPTXM1 

DIFl = PSI { I-l , J+1 )-PSI ( I-l, J-1 )+PSI (1+1, J-1 )-PSI (1+1, J+1 ) 

9 XTERMd.J) = DX*(C0N1*(ZETA(I-1,J)-ZETA(I + 1,J)+ZETA(I-1,J-1)-ZETA 
A( I+l, J-1) )+C0N2*{ V( I , J) *ZETA( I , J) +V( I , J-l ) *ZET A ( I , J-1 ) ) - (U( I , J) * 
BDIF1+U( I, J-1 )*(PSI ( I-l, J)-PSI (I-l, J-2)+PSI ( I+l, J-2S-PSI ( I+l, J) ) )/ 
CC0N3 + CONA*( V( I , J) *(PSI ( I , J+ 1 ) -2 . *PS I ( I , J ) +PS I ( I , J- 1 )) +V ( I , J-1 )* 
D(PS1( I,J )-2.*PSI ( I, J-D+PSI ( I, J-2) ) ) ) 

10 XTERM( I,2)=XTERM( I ,3) 

DO 20 J=2,NPTXM1 
DO 19 I=2,NPTYM2 

DIFl = PSI(I-l.J+l)-PSI(I-l,J-l)+PSI{I+l,J-l)-PSI(I+l,J+l) 

19 YTERM(I,J) = DY*(C0N5*(ZETA(1,J+1)-ZETA(I, J-1)+ZETA(T+1,J+1)-ZETA 
A( I + l, J-1) )-C0N2*(U( I, J)*ZETA( I, J)+U( I +l , J ) *ZETA ( I + 1 , J ) )- (U( I, J )* 
B(PSI( I-l, J)-2.*PSI ( I, J) + PSI ( I + l, J) )+U(I+l, J) *(PSI ( I, J )-2.*PSI ( I +1, 
CJ )+PSI ( 1+2, J) ) )/C0N6 + C0N7*( V( I,J)*DIF1+V(I+1,J)*(PSI(I,J+1)-PSI 
D( I, J-D+PSI ( 1+2, J-D-PSI ( 1+2, J+1) )) ) 

20 YTERM(NPTYM1, J)=YTERM(NPTYM2, J) 

C 
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COMMENT — VERTICAL BOUNDARY WALL POINTS. 

C 

CON3 = CON3/2. 

CONS = CON5*2. 

CON7 = 2.*CON7 
DO 30 I=2.NPTYMi 

XTERMI I, 1) = DX* ( CONI* (ZETAI 1-1,1) -ZETA I !♦ 1,1 H-ZETA( 1-1,2 )-ZET A 
A (1+1,2)) +C0N2*( ZETAd ,1)*V( I, 1) + ZETA( 1,2 )*Vt 1,2) ) - (U(I,D* 

B (PSI ( I-1,2)-PSI { I-l,l)+PSn I+1,1)-PSI ( 1+1,2) )+U( I,2)*(PSI ( I-l, 3) 

C -PSI(I-l,2)+PSI(I + l,2)-PSI(l + l,3)))/C0N3 + C0N4* ( V ( I , 1 ) * ( PS I ( I , 3 ) 
D -2.»PSI( I,2)+PSI ( I ,1) )+V( I,2)*(PSI ( 1 ,4)-2.*PSI(I,3)+PSI( 1,2) ) ) ) 

30 XTERM( I.NOPTX) = OX* ( CONI* ( ZETA II -1 , NOPTX) -ZET A(I +1 , NOPTX ) +ZET A 
A (I-1,NPTXM1)-ZETA(I+1,NPTXM1)) +C0N2*(ZETA{ I,NOPTX)*V( I.NOPTX) + 

B ZETA( I,NPTXM1)*V( I.NPTXMl) ) - ( U 1 1 , NOPTX ) * ( PS I ( I -1, NOPTX ) -PS I { I -I , 
C NPTXMl ) + PSI ( I + 1,NPTXM1)-PSI ( I + 1,N0PTX) ) +U ( I , NPTXMl ) * ( PS I (1-1, 

D NPTXMl )-PS I ( I-1,NPTXM2)+PSI ( I+1,NPTXM2 )-PSI ( I+1,NPTXM1 ) ) )/C0N3 
E + C0N4*( V( I .NOPTX) *( PSI ( I . NOPTX ) -2 . *PS I ( I , NPTXMl ) +PS I ( I.NPTXM2 ) ) 

F + V( I .NPTXMl)* (PSI ( I ,NPTXM1)-2.*PSI ( I , NPT XM2 ) +PS I ( I .NPTXM3) ) ) ) 

DO 40 I=2,NPTYM2 

YTERM( 1,1) = DY*(C0N5*(ZETA( I ,2)-ZETA(I ,1)+ZETA(I+1,2)-ZETA(I + 1 ,1) 
A ) - C0N2*(ZETA(I ,1)*U( I, D+ZETAI I+l, 1) *U( 1+1,1) ) - (U(I,1)*(PSI 
B ( I-l, 1 ) -2. *PS I ( I, 1)+PSI( 1+1,1) )+U( 1+1,1)* (PSI ( I,1)-2.*PSI( 1+1,1)+ 
C PSK 1+2,1) ) )/C0N6 + C0N7*( V( I, 1 )*(PSI ( I-l ,2)-PSI ( I-l.D + PS I ( 1+ 1, 

0 1) -PSK 1+1,2) )+V( I+l ,1)*(PSI (I ,2)-PSI( 1,1 )+PSI ( I+2,1)-PSI ( 1+2, 2) ) 
E ) ) 

40 YTERM( I.NOPTX) = DY* ( CONS* ( ZETA ( I , NOP TX ) -Z ETA ( I , NP TXM 1 ) +ZE TA ( I + 1, 

A NOPTX ) -ZET A ( I + l, NPTXMl) )- COM2* (ZETA (I ,NDPTX)*U(I .NOPTX) +ZETA( I + l 
8 ,NOPTX)*U( I + l, NOPTX) ) - ( U ( I .NOPTX ) * ( PS I ( I- 1 .NOPTX ) -2 . *PS I ( I , 

C NOPTX)+PSI ( I + l, NOPTX) ) +U ( I + 1 .NOPTX ) * ( PS I ( I , NOPTX ) -2 . *P S I (I+l, 

D NOPTX)+PS! ( 1+2, NOPTX) ) )/C0N6 + C0N7* ( V ( I , NOPTX ) * ( PS I ( I -1 , NOPTX )- 
£ PSK I-l, NPTXMl )+PS I ( I+l, NPTXMl )-PS I ( I +1 , NOPTX )) +V ( I +1 , NOPTX )* ( PSI 
F ( I .NOPTX) -PSI ( I .NPTXMl )+PSI ( 1 + 2 , NPTXMl ) -PS I ( I +2, NOPTX) ) ) ) 
YTERM(NPTYMl.l) = YTERM ( NPT YM2 , 1 ) 

YTERM(NPTYM1, NOPTX) = YTERM ( NPT YM2 , NOPTX ) 

C 

COMMENT — HORIZONTAL BOUNDARY WALL POINTS. 

C 

IF(PR08N0 .NE. 0) GO TO 49 
ZETAd, 1) =-2./DY/LAMBDA 
ZETA( I.NOPTX) = ZETAd, 1) 

49 CONI = 2.*C0N1 
CONS = CON5/2. 

DO SO J=3, NPTXMl 

XTERMI 1,J) = DX*(C0N1*(ZETA(1,J)-ZETA(2,J)+ZETA(1, J-1)-ZETA(2, J-1) 
A )+C0N2*(ZETA(l,J)*V(l,J)+ZETA(i,J-l)*V(l, J-1) ) - ( U ( 1 , J ) * ( PS K 1, 

B J + 1 )-PSI ( 1, J-1 )+PSK2, J-D-PSI (2, J + 1 ) ) +U( 1, J-1 )*( PSI ( 1, J )-PSK 1, 

C J-2)+PSI (2, J-2)-PSI (2, J) ) )/C0N3 + C0N4 • ( V ( 1 , J ) * ( PS K 1 , J+ 1 ) -2. * PS I 
0 ( 1, J)+PSK l.J-1) )+V( 1, J-1)*(PSI (1, J)-2.*PSI (1. J-1 )+PSKl, J-2) ) ) ) 

50 XTERMINOPTY, J) = DX«(C0N1*(ZETA(NPTYM1, J)-ZETA(NOPTY, J)+ZETA{ 

A NPTYMl, J-1 ) -ZETA (NOPTY, J-1) ) +C0N2* ( ZET A ( NOPTY , J ) * V ( NOPTY , J ) +ZE TA 
B (NOPTY, J-1)*V(N0PTY, J-1) ) - ( U ( NOPTY , J ) * ( PS I ( NPTYMl , J+1 ) -PS I ( 

C NPTYMl, J-1 )+PS I (NOPTY, J-D-PSI (NOPTY,J+l) ) +U ( NOPTY , J-1 )*( PS K 
D NPTYMl, J ) -PSK NPTYMl, J-2) +PS I (NOPTY, J-2) -PS I (NOPTY, J ) ) )/C0N3+CCN4 
E *( V( NOPTY, J)*( PSK NOPTY, J + 1 )-2.*PSI ( NOPTY , J ) +PS I ( NOPTY , J-1 ) ) +V ( 

F N0PTY,J-1)*(PSI (NOPTY, J)-2.*PSI (NOPTY, J-1 ) + PS I ( NOPTY , J-2 ) ) ) ) 
XTERMd,2 )=XTERM( 1,3) 

XTERM(N0PTY,2)=XTERM(N0PTY,3) 

DO 60 J=2, NPTXMl 

YTERMd, J ) = DY*(C0N5*(ZETA(1, J+1)-ZETA(1, J-1)+ZETA(2, J + 1 )-ZETA (2, 

A J-1 ) )-CQN2* 
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* (ZETA( 1, J)»Ua,J)+ZETA(2. J)*U(2,J) )-(U(l,J)*(PSI(l,J)-2.* 

B PSI(2,J)+PSI(3,J))+U{2,J»*(PSI(2.J)-2.*PSI (3, JI+PSI ( <f, J ) ) )/C0N6 + 
C C0N7*( V( 1, J)*(PSn It J+l )-PSI (1 , J-1 H-PSI (2, J-1 )-PSI ( 2, J + 1 » ) +V(2 t J ) 
D *(PSH2,J + 1)-PSI(2,J-1M-PSI(3, J-D-PSI (3t J+l) ) )) 

60 YTERM(M0PTY, J ) = DY* ( C0N5» ( ZE T A ( NOPTY , J+ U -ZET A ( NOPTY , J-1) + ZE T A 
A (NPTYMl, J+D-ZETAINPTYMI, J-1) ) -C0N2* ( ZETA ( NOPTYt J ) *U ( NOPT Y , J ) + 

6 ZETACNPTYMl, J)*U(NPTYM1, J) ) - ( U I NOPT Y , J ) * ( PS I I NPT YM2 1 J ) -2 , *PS I ( 

C NPTYMlt J)+PSI (NQPTY, J) ) +LI ( NP TYMl , J ) * { P Sit NP T YM3. J ) -2 . *PS I (NPT Y M2 1 
C J)+PSI(NPTYM1,J))) /C0N6+C0N7*( VtNOPTY, J)* tPSI ( NPTYMl , J + l ) -PS I ( 

E NPTYMl, J-1) +PS It NOPT Y, J-1 )-PSI (NOPTY, J + l) )+V (NPTYMl, J ) *( PSI ( 

F NPTYM2,J + 1)-PSI (NPTYM2, J-1) +PS I (NPTYMl, J-1) -PS I (NPTYMl, J + l) ) ) ) 

C 

COMMENT — CORNER POINTS. 

C 

C0N3 = C0N3/2. 

CONS = C0N5*2. 

C0N7 = C0N7+2. 

DIFl = PSI ( 1,2)-PSI (1,1)+PSI (2,1)-PSI (2,2) 

XTERMt 1, I ) = C0N2*(ZETA( 1, 1 )*V( I, 1 )+ZETA(l ,2)*V(1,2))-(U(1,1)»DIF1 
A +U(1,2)*(PSI(1,3)-PSI(1,2)+PS I 

B (2,2)-PSI(2,3))) /C0N3+C0N4*( V( 1 ,1 ) *(PSI(1,3)-2.*PSI(1,2)+PSI(1,1) 
C )+V(l,2)*(PSI(l,4)-2.*PSI(l,3)+PSI(l,2))) 

YTERMt 1, 1 ) = C0N5*( ZETA( 1 ,2 )-ZETA( 1 , 1 ) + ZET A ( 2 , 2 ) -ZE T A ( 2 , 1 ) )-(U( 1,1 
A )*(PSI(1,1)-2.*PSI(2,1)+PSI(3,1))+U(2,1)*(PSI(2,1)-2.*PSI(3,1)+ 

B PSI ( A, 1 ) ) )/C0N6+C0N7* ( V( 1, 1 )*DIF1 
C +V(2,1)*(PSI(2,2)-PSI(2,1)+PSI(3,1)-PSI(3,2))) 

DIFl = PSI ( 1 ,NOPTX)-PSI ( l.NPTXMl )+PSI (2 ,NPTXM1 )-PSI ( 2.N0PTX ) 

XTERMt l,NOPTX) = C0N2* ( ZETA(1,NPTXM1)*V(1,NPTXM1)+ZETA(1,N0PTX)*V 
A ( 1 ,NQPTX) )-(U( l,N0PTX) *01F1+U( L,NPTXM1 )*(PSI(1,NPTXM1)-PSI(1, 

8 NPTXM2 )+PSI ( 2,NPTXM2 )-PSI (2,NPTXM1 ) ) ) / C0N3+C0NA* ( V ( 1 , NOP TX ) * ( P S I 
C ( 1,N0PTX)-2.*PSI ( l.NPTXMl )+P SI ( 1 , NPTXM2 ) ) +V ( 1 , NPTXMl ) * ( PS I (1, 

D NPTXM1)-2.*PSI( 1,NPTXM2)+PSI (1,NPTXM3) ) ) 

YTERMt 1,N0PTX)= CON 5* ( ZE T A ( 1 , NOPTX ) -Z ET A ( 1 , NPTXM 1 ) +Z ET A { 2 , NOPTX )- 
A ZETAt 2, NPTXMl ) ) - ( U ( I , NOPTX ) * ( P S I ( 1 , NOPT X ) -2 . *PS 1 1 2 , NOPTX ) +PSI ( 3, 

B NOPTX))+U(2,NOPTX)*(PSI(2,NOPTX)-2.*PSI(3,NOPTX)+PSI(A,NOPTX)))/ 

C CQN6 + CnN7*( V( 1 , NOP TX ) *D I F 1 + V ( 2 , NOPTX ) * (PSI ( 2 , NOPTX ) -PS I ( 2 , NPTX Ml ) 
D +PSI (3,NPTXM1)-PSI (3.N0PTX) ) ) 

IFtPROBNO .NE. 0) GO TO 70 
ZETAt 1,1)= 0.0 
ZETAt 1, NOPTX) = 0.0 

70 XTERM(1,1)= DX*(C0Nl*(ZETA(l,l)-ZETA(2,l)fZETA(l,2)-ZETA(2,2) ) + 

A XTERMt 1,1) ) 

YTERMt 1,1)= DYMYTERMt l,l)-CON2*(ZETA(l,l)*U(l,l) + ZETA(2, 1)*U(2,1 
A ) ) ) 

XTERMt 1, NOPTX) = DX * ( CON I* t ZETA ( 1 , NOPTX ) -Z ET A t 2 , NOP TX ) +ZE T A ( 1 , 

A NPTXMl)- ZETA(2, NPTXMl) )+XTERM( 1, NOPTX) ) 

YTERMt 1,NDPTX)= D Y* ( YTERM ( 1 , NOPTX ) -C0N2 * t Z ET A ( 1 , NOPTX ) *U ( 1, NOPTX) 

A +Z6TA(2, NQPTX)*U(2, NOPTX) ) ) 

DIF1= PSI ( NPTYMl, 2) -PSI t NPTYMl , 1 ) +PS I (NOPTY, 1 ) -PSI t NOPTY, 2 ) 

XTERMt NOPTY, 1 ) = DX* t CON 1* ( ZE TA t NPT YM 1 , 1 ) -ZE T A ( NOPT Y , 1 ) +Z ET A ( 

A NPTYMl, 2) -ZETA (NOPTY, 2) ) +C0N2* ( ZE TA ( NOPT Y , 1 ) *V ( NOPT Y , 1 ) +Z FT A ( 

B N0PTY,2) »Vt N0PTY.2 ) )-(U(NOPTY, 1 ) *D I F l + U t NOPT Y , 2 ) * ( PS I ( NPTYM1.3 )- 
C PSI (NPTYMl, 2) +PS I (NOPTY, 2) -PSI (NQPTY,3 ) ) ) /C0N3+ CONA* t V ( NOPTY , 1) 

D *(PSI(N0PTY,3)-2.*PSI(N0PTY,2)+PSI (N0PTY,1) )+V( NOPTY, 2)*(PSI( 

E N0PTY,4)-2.*PSI ( N0PTY,3) +PSI ( NOPTY, 2 ) ) ) ) 

YTERM (NOPTY, 1 ) = DY* ( CONS* ( ZETA ( NOPTY , 2 ) -Z ET A ( NOPTY , 1 ) +ZETA (NPTYMl 
A ,2 )-ZETAt NPTYMl, 1 ) )-C0N2* ( ZETA ( NOP TY , 1 ) *J (NOPTY, 1 ) +Z ET A ( NP TYM 1 , 1 ) 
B *U(NPTYM1,1) )-(U(NOPTY,l)*(PSI ( NPTYM2 , 1 ) -2 . *PS I ( NPTYMl , 1 ) +PS I t 
C NOPTY, 1 ) )+U(NPTYMl ,1 )* (PSI (NPTYM3, 1) -2.*PSI ( NPTYM2, 1 ) +PS I ( NPTYMl, 
0 1 ) ) ) /C0N6+C0NT* (V( NOPTY, 1)*I3IFI + V(NPTYM1,1)*(PSI( NPT YM 2,2) -PSI( 
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E NPTYM2» l)+PSI(NPTYMl,l)-PSnNPTYMlt2) ) ) ) 

DIF1= PSI(NPTYM1,N0PTX)-PS1 ( NPTYMl , NPTXMl ) +PS I ( NOPT Y , NPTX M 1) -PS I 
A (NOPTY,f^OPTX) 

XTERM(lMOPTY,NOPTX)= OX* ( CON I * { ZETA ( NPTYM 1 # NOPTX )-ZETA ( NOPTY » NOP TX ) 
A +Z ETA ( NPTYM 1, NPTXM 1)-ZETA< NOPT Y, NPTXMl ) ) + CON2* I ZETA( NOPTY, NOPT X) * 
B VC NOPTY,NOPTX) +ZETA I NOPTY, NPTXMl) *V (NOPTY, NPTXMl) ) -(U( NOPTY, NOPTX 
C )*DIFl^U(NOPTY,NPTXMl)*( PS n NPTYMl, NPTXMl ) -PS I ( NPTYMl , NPTXM2 ) +PS I 
D ( NOPTY, NPTXM2) -PS I (NOPTY, NPTXMl) ) ) /CON3^C 0N4 *( V( NOPTY , NOPTX )* ( PS I 
E (N0PTY,N0PTX)-2.*PSI ( NOP TY , NPTXMl ) +PSI ( NOPT Y , NPTXM2 ) )+V( NOPTY, 

F NPTXMl ) *(PSI (NOPTY, NPTXMl) -2.* PS I ( NOPT Y , N PTXM2 ) +PS I ( NOPT Y , NPTXM3 ) 
G ) ) ) 

YTERM(NOPTY,NOPTX)= 0 Y* (CON5* ( ZETA ( NOPT Y , NOPTX ) -ZET A { NOPTY , NPTX Ml 
A )+ZETA(NPTYMl,NOPTX)-ZETA(NPTYMl,NPTXMl) )-C0N2*CZETA (NOPTY, NOPTX) 
B *U ( NOP TY, NOPTX )-ZETA( NPTYMl, NOPTX ) *U( NPTYMl , NOPTX) )-(U( NOPTY, 

C NOPTX)*( PSI (NPTYM2,NOPTX)-2,*PSI (NPTYMl,NOPTX)+PSI (NOPTY,NOPTX ) ) -f- 
D U( NPTYMl,NOPTX)*(PSI ( NPT YM3, NOPTX ) -2. *PSI ( NPTYM2 , NOPTX ) ^-PS I ( 

E NPTYMl, NOPTX) ) ) /CON6-J-CON7* ( V( NOPTY, NOPTX) *0 1 Fl + V ( NPT YMl , NOPTX ) * 

F (PSI (NPTYM2, NOPTX) -PSI ( NPT YM2 , NPTXMl ) *PS I ( NPTYMl , NPTXM 1 ) -PS I ( 

G NPTYMl, NOPTX) ) ) ) 

OMMENT — INITIALIZE PRESSURE FIELD. 

DO 200 J^I, NOPTX 
DO 200 1=1, NOPTY 
200 P( I , J )= 1.0 
ITKONT= 0 

202 ITKONT= ITKONT+1 
JAIL= .FALSE. 

DO 210 l=I,NPTYM2 
DO 210 J=l, NOPTX 

204 IF (J .GT. 2) GO TO 205 
PX= P( I, J+1)-XTERM( I, J) 

GO TO 207 

205 PX= P( I, J-l)+XTERM( I, J) 

207 PNEW= 0.5*(PX+P(I+1,J)+YT£RM{I,J)) 

IF (JAIL) GO TO 209 

IF (ABS( (PNEW-P(I , J) ) /PNEW) .GT. PCT ) JAIL = .TRUE. 

209 P ( I , J ) = PNEW 

210 CONTINUE 

DO 220 I=NPTYM1, NOPTY 
DO 220 J=l, NOPTX 

214 IF (J oGT. 2) GO TO 215 
PX = P(l,J+l)-XTERM(I,J) 

GO TO 217 

215 PX = P( I, J-l)^XTERMn,J) 

217 PNEW = 0.5*(PX+P( I-l, J)-YTERM(I , J) ) 

IF (JAIL) GO TO 219 

IF (ABS( (PNEW-P(I,J))/PNEW) .GT. PCT) JAIL =.TRUE. 

219 P( I , J ) = PNEW 

220 CONTINUE 

IF ( .NOT. JAIL) GO TO 229 
IF(ITKONT .LT. 500) GO TO 202 
WRITE (6,225) 

225 FORMAT! 1H14X,50H*»***N0 CONVERGED SOLUTION IN 500 I TERATI ONS.** *♦♦ ) 
A) 

229 WRITE (6,230) 

230 FORMAT! 1H18X, 18HPRESSURE FIELD, P*) 

IF (RE .LE. 1.) WRITE (6,231) 

231 FORMAT! IH+26X, IH« ) 

WRITE (6,232) ITKONT 
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232 FORMAK 1H08X.19HN0. QF ITERATIONS =I4/1X) 

DO 235 I=l,NOPTY,YINDX 
235 WRITE (6fl06) { P ( I , J ) , J=1 , NOPTX , X I NDX ) 

106 FORMAT(9X,15F8,4) 

IF (JAIL) GO TO 202 
DO 240 J=1,N0PTX 
DO 240 I=1,N0PTY 

240 PNORM(I,J)= P( I ,J)-P( IREF, JREF) 

WRITE I6t241) 

241 FORMAT! 1H18X,26HN0RMALIZED PRESSURE FIELD, /IX) 

DO 245 I=1,N0PTY, YINOX 

245 WRITE (6,106) ( PNORM( I , J ) , J= 1 , NOPTX , X I NDX ) 

DX03 = DX/3. 

NPTXPl = NOPTX+i 
DO 250 I=1,N0PTX 
II = NPTXPl - I 

IF (PN0RM(1,II) .LT. 0.) GO TO 252 
250 CONTINUE 
252 II = II+l 
LEREG = 0. 

IF (M0D(II,2) .NE. 0) GO TO 255 

LEREG = (PNORM( 1,11 )-i-PNORM(l,m-l) )*DX/2. 

II ^ II+l 

255 IF (II. EQ. NOPTX) GO TO 263 
DO 260 J=II,NPTXM2,2 

260 LEREG = ( PNORM ( 1 , J ) +4 . *PNORM ( 1 , J+ 1 ) +PNORM ( 1 , J +2 ) ) *0X03 + LEREG 
263 WHSURF = 0* 

DO 265 J=1,NPTXM2,2 

265 WHSURF = WHSURF + ( PnORM ( 1 , J ) +4 . * PNORM { 1 , J +1) + PNORM ( 1 , J +2 ) ) *DXO 3 
WRITE (6,505) LEREG, WHSURF 

505 FORMAT! 1HK8X,35HNET FORCE FOR LEADING EDGE REGION =1PE1 3. 5/23X , 15H 
AWHOLE SURFACE =E13.5) 

IF (NOPLOT -EQ. 0) RETURN 
DO 333J=1, NOPTX 

333 CALL BCDUMP(PN0RM(1, J) ,PNORM(NOPTY, J) ,1 ) 

RETURN 

END 
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$IBFTC WFIELO LIST, DECK 

c 

COMMENT — W* VELOCITY PROFILES 
C 

DIMENSION PS I (65,65) ,W( 65,651 , DPSI DY { 65 , 65 ) , OPS I DX ( 65 , 65 ) , 
A WPCT(65,65) 

REAL LAMBDA 
INTEGER XINOX,YINDX 
LOGICAL INDKT 

READ ( 5, 3 ) LAMBDA,RE,OMEGA,PCT,NOPTY,NOPTX,NOL INYtNOLINX 
3 F0RMAT(4F6.0,4I3) 

NPTXMl = NOPTX-1 
NPTYMl ^ NOPTY-1 
XINDX = NPTXM1/(N0LINX-1) 

YINDX = NPTYMl/(NQLINY-l) 

BSODSQ= 1. /LAMBDA** 2 
DX = 1*/FL0AT(NPTXM1) 

DY = 1*/FL0AT(NPTYM1) 

DXSQ= DX*DX 
OYSQ= OY*DY 
TWODX= DX+DX 
TWOOY= DY+OY 

WPF = 2.*(1./DXSQ + BSQDSQ/OYSQ)/RE 
ONMOM = - OMEGA 

OMOWPF= OMEGA/WPF 
C 

COMMENT — SET UP STREAM FUNCTION DISTRIBUTION 
C 

READ (5,5) ((PSI(I,J) , I=1,N0PTY) , J=1,N0PTX) 

5 FORMAT (9F8.0) 

DO 8 J=1.N0PTX 
DO 8 I=1,N0PTY 
8 PSI(I,J)= PS! ( I ,J) *l.E-4 
DO 12 J=2, NPTXMl 
DO 12 1=2, NPTYMl 

OPSIDY( I, J)= (PSI ( I-l, J)-PSI ( I+l, J) )/TWODY 
12 DPSIDX(I,J)= (PSI ( I , J+1 )-PSI ( I t J-1 ) )/TWOOX 
C 



COMMENT — INITIALIZE W* FIELD 
C 

DO 14 J=2*NPTXM1 
W(NOPTY,J) = 0.0 
DO 14 I=1»NPTYM1 
14 W( I ,J)= -1.0 
W(l,l)= -1.0 
W(1,N0PTX) = -1.0 
DO 16 I=2,N0PTY 
W(I,1)= 0.0 

16 W( ItNOPTX) = 0.0 

17 WRITE (6,28) 

28 F0RMAT( IHl) 

READ (5,18) DPDZ 

18 FORMAT! F8.0) 

ITKONT = 0 

C 

COMMENT — ITERATE W* FIELD FOR GIVEN DP/DZ VALUE. 

C 

20 INDKT = .FALSE. 

ITKONT = ITKONT 4 1 
DO 22 I=2,NPTYM1 
DO 22 J=2,NPTXM1 

WF= DPSIDY( I,J )*( W( I , J+1)-Wn ,J-1) ) /TWODX - DPS I DX ( I , J ) *( W ( I-l , J ) 
A-W( I + 1,J) )/TWODY+DPDZ-( ( W ( I , J4-1 ) +W ( I , J- 1 ) ) /DXSQ+BSQDSO*( W ( I-H,J ) + 
BW( I-l, J ) )/OYSQ)/RE 
WNEW = ONMOM*W(I,J) - OMOWPF*WF 
WPCT(I,J)= ( W( I ,J)-WNEW)/WNEW 
IF ( INDKT) GO TO 22 

25 IF( ABS(WPCT(I,J) ) .GT. PCT) INDKT - .TRUE. 

22 W( I,J)= WNEW 

IF (INDKT) GO TO 20 
C 

COMMENT — CALCULATE Q,NET FROM CONVERGED W* FIELD. 

C 

QNET = 0. 

DO 60 I=2,NPTYM1,2 
DO 60 J=2,NPTXM1,2 

60 QNET = QNET+16.*W(I,J)+4.*CWCI-1,J)+W(I+1, J)^W(I,J-1)+W(I,J+1)) 

A +wn-l, J-l)+W( I-l, J+l)+W( I + I,J-1)^“W( I + 1,J + 1 ) 

QNET = QNET*0X*DY/9. 

C 

COMMENT — OUTPUT 
C 

WRITE (6,26) LAMBDA, RE, OMEGA, PCT, ITKONT, NPTXM1,NPTYM1, DPDZ, QNET 

26 FORMAT! 9X,8HLAM80A =F6. 3 , 5X, 4HRE =F7. 1 , 5X , 1 9HREL AXAT ION FACTOR 
A=F5.2/9X,24HMAXIMUM RELATIVE ERROR =F6. 3 , 5 X, 22HNUMBER OF ITERATION 
BS =I5/9X,37HNUMBER OF INCREMENTS IN X-DIRECTION =I4,5X,16HIN Y-DIR 
CECTION ^I4/9X,9HDP*/0Z* =F7. 3 , 6X, 6HQ* , Z = 1PE13.5/1H08X,20HW* VELO 
DCITY PROFILE, /IX) 

DO 30 I==1,N0PTY, YINDX 
30 WRITE (6,33) ( W ( I , J ) • J= 1 , NOPTX, XINDX ) 

33 FORMAT! 1X,15F8.4) 

GO TO 17 
END 
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SIBFTC GRAPH LIST, DECK 

COMMON LAMBDA,NOPTY,NOPTX,NOYPLSf NOXPLS,DY,DX,PNORM 
COMMON /SKALF /XSCALF , YSC ALF , Z SCALF 
DIMENSION Z( 13000) ,PNORM(65, 65) 

EXTERNAL PSURF 
REAL LAMBDA 

READ (5,4) L AMBDA, NOPTY, NOPTX, NO YPLS, NOXPL S, XSCALF, YSC ALF, PSCALF 
4 FORMAT(F6,0,4I3,3F6.0) 

ZSCALF = PSCALF 

DY= l./FLOAT(NOPTY-l) 

DX= l./FLOAT(NOPTX^l) 

NOPTS = 3*(N0PTX*N0YPLS+N0PTY*N0XPLS) 

DO 10 J=l, NOPTX 

10 CALL BCREAD(PNORM(1,J),PNORM(NOPTY, J) ) 

CALL PL0T3D ( 0 1 0. , LAMBDA, Z , NOPTX, NOYPL S , NOXPLS , NOPTY, PSURF , 

A .TRUE.) 

CALL ROTATE ( 0. ,0. , 45 . , . F ALSE . ) 

CALL ROTATE (0. ,35. ,0. , .TRUE. ) 

STOP 

END 


SIBFTC SURFAC LIST, DECK 
FUNCTION PSURF(I,J) 

COMMON L AMBDA, NOPTY, NOPTX, NOYPL S, NOXPLS, DY ,DX,PNOR M( 65, 65 ) 
I 1= NOPTY + 1 -I 
PSURF = PNORMlII,J) 

RETURN 

END 


RMAP 

• TCRE A 




rNTRY 

BCRFAT 


AO 

SA VF 

1 



Cl. A 

T ,4 

GET FIRST ARG. 


LTQ 


GET SECOND ARG. 


TLQ 

*4-? 

CO'^PARE 


XC A 


IF 2ND LESS EXCHANGE 


STO 

M r> 

STORE SMALLEST ARG 


ADD 

S YSONF 

ADD 1 


STA 

STH 

STORF FOP MOVE 


S'JR 

TFMP 

COMPIJTF COUNT 


STA 

1 XI 

STORE FOR MOVF 


A YC 

UNC-T ,4 

LOCATE UN^S LIKE F IV CALL 


SXA 

SYSL0C,4 

AND S AV F IN SYSLOC 


CALL 

. . 

SET IJP RFAD 

Q- AT 

TSX 

. . F I nc\ 4 

READ RECORD 


TSX 

• . FTCK,4 

CHECK READ 

I 

A XT 

, 1 

PICK UP COUNT L EFT 


TXL 

LASTCa 

IS ONLY I RFC L EFT 

T 

A XT 

,4 

RFC CNT 


CL A 

• . 5P0B4-? ,4 

MOVE WHROS 


STH 


TO STORE 


TI X 

,l ,1 

DECR. COUNT 


TXI 

^4-1 ,4,-1 

DFCR. RFC COUNT 


41 



C I 4 

TXM 

STO-1 ,4 


SXA 

I Yl ,1 


A 

D E AO 

M st: 

TXL 

PONE ,1 


A XT 

0QNF,4 


SXA 

L ASTC-1 t4 


SCO 

CKI P4,l 


TR A 

r X4 

O'^MF 

AXT 

RF ao,4 


SXA 

1. ASTC-l T V 


AX^ 

”? P r4 



CKI R4,4 


R F TIJP M 

RCPE AO 

UNC 

P7E 

. UN"'E. 

T" M p 

P7P 



*Lnf R 
rr.^n 


CP, P=C CUUNT 
NO SAVE REMAINING COUNT 

on reao next recur n 

ANV MORP WORDS 

YFS STORE TP FX IT NEXT 

TIME 

SET RFC CNT = NO WORDS LEFT 
GO PROCESS RFCGRO 
RESTORE EX IT 

RESTORE REC CNT 


Aon or UNIT s 


\ 


f f B M A ^ 

. BC RWO 



c NITRV 

. . OP OB 


r ^}TR Y 

. . 3CPO 


NTPY 

. . BCWO 

SI 7<“ 

S^T 

P? 

. . « C P 

C LA 

POO 


TO A 

- 4- P 


CAL 

PTO 


SXA 

LK. OK, 4 


CALL 

. , PIPSf SFL) 


ORG 


s: L 

T OPT 

. . BPOB, ,5I7F 


L XA 

LK. OP ,4 


TR A 

1 ,4 

PON 

PON 

t »o 

PT^ 

p TM 


. , RROO 

BSS 

SI 7.F 

LX. OP 

LOT R 



R on 


RFCORD SIZE 
READ ENTRY FOR BCR FAD 
GET CORRECT AR G FOR ,EIOS 
wRiTF entry for BCDUMP 
SAVE IP4 

$PT UP READ OR WRITE 

I/O COMMAND 
RESTORF 4 


I/O BUFFER 


42 



APPENDIX E 


SAMPLE PROBLEM 

A sample problem, which was solved on the Lewis IBM 7044-7094 direct-couple 
system, is included to show the user how to set up the input data cards and to show what 
output can be expected. For this problem LAMBDA = 1 (square groove), RE = 100, 
OMEGA = 1.3, PCT = 0. 001, NOPTY and NOPTX = 29, and NOLINY and NOLINX = 15. 
The reference mesh point for the normalized pressure field calculations is (29, 15), and 
the punched output for the w*- velocity calculations and normalized pressure field plot is 
required. The control words PROBNO and FLUFLG are zero in this problem. Results 
for this case are discussed in reference 1. The input data card for the ISGPSM program 
should be punched as 


1 8 

9 16 

17 

24 

25 

32 

35 

38 

41 

44 

47 

50 

53 

56 

59 

62 

63 

'72 

73 

80 

1.000 

100.0 


1. 30 


0.001 

29 

29 

15 

15 

.29 

. 15 

1 

1 

0 

0 






The printed output for this problem is shown in table I. The number of iterations 
printed in table 1(a) is the number of iterations required for convergence of the stream 
function to within the designated relative change. It should be noted that the values of the 

4 

stream function have been multiplied by 10 before printout to facilitate writing the out- 
put format of the stream function and to present more s^nificant figures in the individual 
values. In table 1(b) the values of the vorticity at the four corners are those calculated 
from the equations for the vertical stationary walls. If the values at the four corners 
based on the equations for horizontal walls are desired, they are = 0 for the lower 
stationary wall and = -2/XAy* for the upper moving wall. The static pressure field 
is printed in table 1(e) along with the number of iterations to convergence to within the 
designated relative change. The relative change for convergence of the pressure field is 
equal to that of the stream function. The normalizing factor for this particular problem 
is 1. 0025 and is the entry found in the eighth column of the bottom row. This number is 
used to calculate the normalized pressure field distribution (table 1(f)). Also shown in 
table 1(f) is the leading edge force and net force actii^ on the moving surface. 

Using the punched output of the stream function from the ISGPSM program, the w*- 
field and Q* can now be calculated using the WFIELD program. The additional data 

Zi 
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TABLE r. - SAMPLE PROBLEM OUTPUT 


(a) Stream function 
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(b) Vorticity 
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TABLE I. - Concluded. SAMPLE PROBLEM OUTPUT 
(d) V*- Velocity profile 
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required are that LAMBDA = 1, RE = 100, PCT = 1. 30, NOPTY = 29, NOPTX = 29, 
NOLINY = 15, NOLINX = 15, and 3P*/9z* = -0. 115. The first WFIELD program input 
data card is punched as 
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73 

80 

1.000 


100.0 
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29 
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15 


15 







This card is followed by the deck of cards containing the stream function values of 
ISGPSM and the card with the 9P*/9z* value. 


1 8 

9 

72 

73 

80 

- 0. 115 






K more 9P*/9z* values are to follow for a given problem, they must be punched 
one to a card according to this last format. 

Table II contains the one page of output per 9P*/9z* value from the WFIELD pro- 
gram. The number of iterations printed is the number of iterations of the w*- velocity 
profile from initial value to converged value within the prescribed maximum relative 
error. The net volume flow Q* associated with the 9P*/9z* value of -0. 115 in this 
problem is 0.0003128 and is printed on the same line as the 9P*/9z* value. 

The computer execution time on the Lewis Computer for this sample problem was 
approximately a half a minute for each program. 
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The three-dimensional plot of the normalized pressure distribution is obtained from 
program GRAPH using the column binary punched output of the pressure distribution 
from the ESGPSM program. The additional information required is that LAMBDA = 1; 
NOPTY = 29, NOPTX = 29, NOYPLS = 15, NOXPLS = 15, and XCALF, YSCALF, 
PSCALF each equal 6. 0. The first program input data card is punched as 


7 9 

10 

12 

13 

15 

16 

CO 

19 

24 

25 

30 

31 

36 

37 

72 

73 

80 

29 


29 


15 


15 


6.00 


6.00 


6.00 






After this is the deck of cards containing the pressure distribution values. 

The output for this particular problem is the graph found in figure 15. The maximum 
value of P* = 1. 2501 is found at the corner where the movir^ wall meets the leading 
edge and the minimum value of P* = -0. 3304 is near the corner formed by the moving 
wall and the trailing edge. Note also the shallow pressure drop or relative minimum in 
the vicinity of the vortex center. 
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Figure 1. - Spiral groove pumping sea) model for limiting case of Figure 2. - Streamlines In groove cross flow plane, 

land clearance. (+ designates vortex center. ) 



Figure 3. - Mesh point representation of x -y plane Figure 4. - MAIN program, 

flow field. 
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Figure 5. - Stream function and vorticity calculations. 
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Figure 6. - Stream function initial distribution. 
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Figure 7, - Vorticity initial jistribution. 
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Figure 10. - Concluded. 
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Figure 12. - Three-dimension pressure plot. 






Figure 13. - Finite clifference approximation of f.c iDCL-nilary points. 



Figure 14. - Overlay structure of spiral grooved pumping sea! 
model programs. 
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